Many body localization in Heisenberg XXZ magnet in a random field 
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We numerically investigate Heisenberg XXZ spin- 1/2 chain in a spatially random static magnetic 
field. We find that tDMRG simulations of time evolution can be performed efficiently, namely the 
dimension of matrices needed to efficiently represent the time-evolution increases linearly with time 
and entanglement entropies for typical chain bipartitions increase logarithmically. As a result, we 
show that for large enough random fields infinite temperature spin-spin correlation function displays 
exponential localization in space indicating insulating behavior of the model. 
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I. INTRODUCTION 

Interesting open problems in solid state theory are fre- 
quently connected to strongly correlated systems. While 
the physics of noninteracting particles is well under- 
stood much less in known about many body interact- 
ing systems. Analytical solutions are typically not pos- 
sible and numerical calculations are notoriously difficult 
precisely due to strong correlations. In the theory of 
quantum information quantum correlations, or entangle- 
ment, are one of the central objects of study. Recently, 
ideas of quantum information inspired a new numerical 
method, called time-dependent density-matrix renormal- 
ization group (tDMRG)ii^. The method, originally pro- 
posed for simulation of time evolution of quantum sys- 
tems on classical computers, has been shown rigorously 
to be efficient in the number of particles at a fixed evolu- 
tion time^. But the more relevant question is what is its 
time efficiency, i.e. how the complexity grows with sim- 
ulation time? Numerical experiments showed that the 
method is efficient only in some rather special cases^, 
while in general it fails^ due to fast growth of entangle- 
ment with time. An open question is whether tDMRG 
can nevertheless be used efficiently for some nontrivial 
interacting many particle system and for generic initial 
conditions? Our aim is to answer this question in af- 
firmative. We show that tDMRG is time efficient for an 
interacting Id Heisenberg chain in a disordered (spatially 
random) magnetic field. This efficiency allows us to cal- 
culate infinite temperature correlation functions for large 
chains of the order of a hundred sites and showing that 
the many-body states in disordered Heisenberg model are 
localized, at least at large enough random fields. 

While a single particle Anderson localization is well 
understood, for a review see^, the interplay between dis- 
order and (strong) interactions in the onset of localiza- 
tion, as manifested e.g. in vanishing d.c. transport coeffi- 
cients, is a subject of an ongoing research, see e. g. refsji 
and references therein. The simplest interacting situation 
is that of two particles. It has been studied for the first 
time in* and shown that the localization length can dras- 
tically increase in the presence of interaction. This has 
been confirmed by many subsequent studies. Situation 



for many interacting particles, e.g. for finite densities, is 
much less clear. It has been explored ii*^ by calculating 
time evolution of wavepackets, showing that the center of 
mass extension of the wavepacket grows logarithmically 
with time. Later, the influence of the disorder on the en- 
tanglement has been studied for single particle states— 
and for quantum computer simulating single particle lo- 
calization^^. Disordered Heisenberg model and entangle- 
ment properties of its eigenstates has been studied in^-^ 
where a transition in level spacing distribution from Pois- 
sonian for no disorder, to Wigner-Dyson distribution of 
random matrix theory, and back to Poissonian in the 
case of localization, has been observed as the disorder 
amplitude is increased. Spectral statistics for interacting 
disordered system has been also studied iv^, sugesting 
the existence of localization for sufficiently strong dis- 
order. Characterization of metal-insulator transition in 
disordered systems (in absence of interaction) in terms of 
spectral statistics have been studied in numerous works, 
see, e. g. reffli and references therein. Similar tran- 
sitions in spectral statistics have been observed also in 
interacting systems, see e.g. refs.— . Ini^ a relation be- 
tween generalized entanglement and inverse participation 
ratio has been obtained for eigenstates of a disordered 
Heisenberg model whereas Meyer- Wallach entangle- 
ment has been calculated for random states localized on 
M random or adjacent computational states. Localiza- 
tion in many body system can be obtained also by con- 
structing special on-site disorder—. 



II. NUMERICAL METHOD 



We will use tDMRG method to calculate time evo- 
lution of pure states^, \i^{t)) = U{t)\^[Q)), or time 
evolution of operators^, 0{t) = W{t)OU{t). Matrix 
product states (MPS) are used to represent pure states, 
W = ti'C^r • • • ■ • Sn), where |si . . . s„) are 

computational basis states with each sj taking two values 
Sj € {0, 1}, i.e. local dimension is d = 2. For operators 
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a matrix product operator (MPO) is expansion 

where a basis of Pauli matrices is used for each site, 
that is each Sj can now take four different values, Sj € 
{0,x, y,z} (local dimension is d = 4). The advantage 
of MPS/MPO representation is that the transformation 
acting on the neighboring spins can be done locally, that 
is by transforming only two adjacent matrices. Short 
time propagator U{t) = exp {—iHt) generated by near- 
est neighbor Hamiltonian H is decomposed using a 2nd 
order Trotter-Suzuki formula into a series of one and 
two qubit (qudit) operations. After each application of 
a two qubit gate the dimension D of the two matrices 
involved increases by a factor of d. In order to prevent 
the growth of matrix size with time one truncates their 
size using a singular values decomposition, keeping only 
the largest singular values. Truncation after application 
of a single two qubit gate Uk introduces truncation er- 
ror equal to the sum of squares of the discarded sin- 
gular values, 7j{Uk) = 1 - J2f=i f^'ji^k), if A'j (f^fe) are 
decreasingly ordered singular values {i.e. Schmidt co- 
efficients) of the bipartition with the cut being on the 
bond affected by the gate Uk. Total truncation error 
after application of a series of gates, U{t) — J|j, Uk, is 
then the sum of individual errors r]tot{t) = X^fc'K^fc)- 
Truncation error r]tot{t) scales with the time step r as 
'7tot(0 By using Trotter-Suzuki factorization of U{t) 
we also introduce Trotter error. For our choice of 2nd or- 
der formula the error in fidehty 1 - KV'tDMRcl'/'oxact)!^ 
scales as esc T^(i/T)^ — tH"^ (note that the "phase" er- 
ror KV'tDMRolV'oxact) " 1| scalcs as (X T^i). If one starts 
evolution from a product state/operator at i = 0, which 
is always the case in our simulations, the error is ini- 
tially dominated by the Trotter error but for larger times 
the precision of tDMRG is eventually determined by the 
truncation error ?7tot- In the following we are going to 
focus solely on the truncation error. 

Throughout this paper, when calculating evolution of 
pure states, we start with a random product state, i.e. 
IV'(O)) = - • •(8' IV'n), where \tpj) is a random state of 

j-th qubit corresponding to a random point on its Bloch 
sphere. Note that in numerical simulations, in order 
to obtain an infinite temperature behaviour, we average 
over initial random states. When simulating Heisenberg 
evolution of operators - which will be used for computa- 
tion of spin-spin correlations - the initial operator will be 
chosen to be the spin projection at a quarter of the lattice 
0(0) = s'^n/i- '^^^ initial state is therefore in both cases, 
of operator or pure state dynamics, separable. However, 
time evolution is expected to produce entanglement. In 
all numerical computations, we also average over random 
realizations of disorder. 

The system we study is a Id spin-1/2 Heisenberg XXZ 



model in a random magnetic field, 

n — 1 n 

^ = E (^^^^+1 + 44+1 + + E h,s). (2) 

where s" are canonical spin 1/2 variables. The mag- 
netic field will be chosen randomly and uniformly in the 
interval [—h,h]. The case of A = is a special case 
and should be clearly distinguished from nonzero A. Us- 
ing Wigner-Jordan transformation, Hamiltonian ^ with 
A = can be transformed to a bilinear fermionic system 
X]j(Cj+iCj + h.c. + hjUi) with rii = c\ci which represents 
the model of noninteracting spinless fermions with diag- 
onal disorder known to exhibit Anderson localization in 
Id. The case A ^ introduces the interaction or cor- 
relations among electrons (through Arij/ij+i) which can 
qualitatively change properties of the system. The aim 
of this letter is twofold: (i) to study localization in an 
interacting disordered system, and (ii) to show that time 
evolution with tDMRG, as opposed to non-integrable and 
non-disordered situation, is efficient for such a system. 

The computational complexity of simulating quantum 
evolution on a classical computer using tDMRG is de- 
termined by the growth of bipartite entanglement. If 
bipartite entanglement increases with time we have to 
increase dimension of matrices D with time in order to 
prevent truncation error 77tot(i) from growing. Because 
the number of computer operations in tDMRG scales as 
^ it is crucial to know how fast do we have to increase 
Dl We are going to study necessary dimension Ds{t) in 
order for the truncation error at time t to be less than 
e. If the necessary D^{t) grows linearly with time we 
say that the simulation is efficient and if it grows expo- 
nentially simulation is inefRcient. For a quantum chaotic 
many body system D^{t) grows exponentially with time^. 
Before looking at (t) let us have a look at bipartite en- 
tanglement. Note that the z-component of the total spin 
= Sj is a conserved quantity [H, S^] = 0. In the 
following we will consider T — > oo limit of the model, so 
that random states will on the average represent states 
with zero polarization, = 0. In the fcrmion repre- 
sentation this corresponds to the half-filled band with 
n/2 spinless particles where the effect of interactions is 
strongest. In the other extreme case S*^ — n/2 — 1 the in- 
teraction A term is constant and therefore the problem is 
equivalent to Id Anderson model of localization without 
interaction. 



III. RESULTS 
A. Entanglement entropy 

The amount of bipartite entanglement can be mea- 
sured by von Neumann entropy Si{t) = tr (p^i log2 pa), 
or in general by the Renyi entropy, of the reduced density 
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FIG. 1: Entanglement entropies So.5{t) and Si{t) for pure 
state evolution (a), and operator space evolution (b) [see text 
for details]. We compare cases of A = 0.5 (growing curves) 
and A = (saturating curves), for n = 50, average over 
100 (a), 21 (b, for A = 0.5), 1000 (b, for A = 0) disorder 
realizations with magnetic field magnitude h — 5. In insets 
we compare logarithmic growth of Si{t) in disordered case 
(full curves, same data as in main frames) to linear growth 
in the case of staggered magnetic field hj — ( — 1)-'5/\/3 and 
A = 0.5 (dashed curves). 

matrix 

^aW = i^^^^, ti-B\m){m\- (3) 

1 — a 

We shall always consider bipartite cut with the first m, 
or last n — m, qubits constituing subspace A, or B, re- 
spectively. In the case of operators, expansion coeffi- 
cients in the basis of Pauli matrices ([T]) are treated as 
expansion coefhcients of a super-ket in a Hilbert space 
of dimension 4", for which the operator- space entangle- 
ment entropy is then calculated. For the discussion of 
approximability of states with MPS form and its rela- 
tion to Renyi entropies sesi^. We show in Fig. [T] time 
dependence of Renyi 6*0.5 (i) and von Neumann Si{t) en- 
tropies. For simulation of pure states, Fig. [1^, the en- 
tropies can be seen to grow logarithmically with time in 
the interacting case A = 0.5, whereas they saturate to 
a constant for A = which seems consistent with an 
Anderson quasi-particle localization. Entropies in Fig. [T] 



are shown for the 'worst case' bipartition only, i.e. for 
m which maximizes them, whereas results are qualita- 
tively equivalent for half- half bipartition (m — n/2), or 
average bi-partition (average over m = 1 . . . n — 1). For 
evolution of operators, Fig. [T}d, in the interacting case 
A = 0.5, 5*0.5 (t) again grows logarithmically with time 
whereas the growth of Si{t) is slower than logarithmic, 
perhaps saturating. In fact, saturation of operator space 
entropies Sa (t) , for large enough a > a* , seems a proba- 
ble interpretation of our results since it is compatible with 
localization of correlation functions reported below (see 
Fig. [5]). Slower growth of entanglement entropies (or even 
saturation) for operators as compared to faster growth for 
pure states is consistent with a qualitatively similar find- 
ing for a homogeneous transverse Ising model^ . In insets 
of Fig.[l]we compare our data to tDMRG simulation for a 
non-disordered model ([2]) in a staggered field of compara- 
ble strength hj = {—iyh/y/3, and there we find a linear 
growth of Sa(t) (X t, for pure states and operators, which 
is consistent with an exponential increase of (t) found 
for non-integrable models^. Note that there is in general 
no simple relation between the behavior of pure state 
entanglement entropy and operator-space entanglement 
entropy. Only in the special case of rank-one projection 
operators O = \^p){'lp\ we find a simple relation, namely 
that operator-space von Neumann entanglement entropy 
of equals two times entanglement entropy of \ip). 

Hence the operator-space entanglement entropy is also 
not equivalent to an entanglement of a mixed state. 

As we have seen, entropies grow at most logarithmi- 
cally for a disordered field. This gives us hope that the 
evolution with tDMRG is in fact efficient, meaning that 
D^(t) grows polynomially with time. This is indeed the 
case as shown in Fig. [2] The necessary dimension of ma- 
trices grows linearly with time therefore the simulation 
of disordered Heisenberg chain is efficient, both for pure 
states and for operators. This must be contrasted to the 
other known efficients case of tDMRG, namely integrable 
transverse Ising model^ , where the simulation is efficient 
only for operators. 

We note that in the numerical results presented above 
there is no significant finite size effects. In order to 
demonstrate that, we plot in Fig. [3] the growth of Renyi 
entropy >S'o.5(t) for different system sizes n, and observe 
significant finite size effects only for size smaller than 30. 

B. Infinite temperature correlation function 

After establishing that tDMRG simulation of a disor- 
dered Heisenberg chain is efficient we want to calculate 
some physically relevant quantity. We choose an infinite 
temperature spin-spin correlation function C(r, <), 

C^('^,0 = 2-"tr{4/4(t)</4+,}, (4) 

which is computed from MPO representation of Heisen- 
berg dynamics ■3^^/4(0 • By calculating C(r,t) we can di- 
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FIG. 2: Dimension -De(t) for pure state (top) and operator 
evolution (bottom). Data for various e clearly indicate linear 
growth with time. A = 0.5, n = 50, and a single disorder re- 
alization with h = 5 (the same realization for all data points) . 
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FIG. 3: Renyi entanglement entropy growth 5'o.5(t) for pure 
state evolution and different system sizes n. We show nu- 
merical results, averaged over 100 realizations of disorder, for 
n = 6, 10, 30 and 50, with short-dashed, dotted, long-dashed, 
and full curves, respectively. Note that n = 30 and n = 50 
cases almost overlap. In all cases A = 0.5 and h = 5 and we 
show the average entropy over all bipartite cuts in order to 
have smoother data. 



rectly assess the many-body localization. In particular, 



C{r,t — > oo) gives direct information on d.c. transport 
properties, i.e., whether the d.c. spin diffusion constant 
is finite, or equivalently for fermions, whether the sys- 
tem is a conductor (normal resistor) or Anderson insula- 
tor (at finite T). First, we are going to consider non- 
interacting system A = 0. The reason to study this 
case is to verify recently obtained bounds on the prop- 
agation of correlations in disordered systems^^. It has 
been proven that the correlations are exponentially sup- 
pressed outside an effective logarithmic lightcone. Note 
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FIG. 4: (Color online) Spin-spin correlation function C(r, t) 
(logj^o color code) for non-interacting case, A = 0, n = 500, 
average over 1000 disorder realizations with h — 1. Solid 
curves are isolines at 10"*, 10"^, 10"" and 10"^, from right 
to left. 

that since this case A = is equivalent to noninteracting 
fermions the necessary dimension of matrices D actually 
saturates^ at Dg{t) = 4, meaning that the calculation is 
very efficient and chains of thousands of spins are eas- 
ily accessible. From the calculated correlation function 
shown in Fig.[¥]we can see that the correlation function 
decays exponentially in space, and is frozen in time after 
sufficiently long time, |C(r, i)| < K exp{—r/£^), for some 
K,(^ > 0. This means that the rigorous estimate^*' is in 
fact over-pessimistic and can perhaps be improved^!. We 
actually observe three regimes: (i) for small times cor- 
relations propagate balistically, seen as a linear growth 
of isocurves C{r{t),t) — const, (ii) After that we have 
a logarithmic propagation of correlations, reflected in a 
log-shape of isocurves. (iii) For large time, localization 
sets in and the spatial correlation function gets its asymp- 
totic exponential shape. The crossover times between the 
regimes increase with decreasing value of correlation on 
the isocurve. 

In Fig. [5] correlation function is shown for the inter- 
acting system, A = 0.5. Interestingly, we again observe 
localization with qualitatively similar structure of corre- 
lation function C{r, t) as in the case A = 0. Note that 
numerical simulation is in this case much harder. For the 
case shown the dimension of matrices has been D = 64 
which, as can be inferred from Fig. [21 means truncation 
error ?7tot(^ ~ 40) « 10"'^. Because of truncation errors, 
for large times t > 20 and distances r > 10 one can 
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FIG. 5: (Color online) Correlation function C{r,t) for n — 50, 
A = 0.5, averaged over 21 disorder realization with h — 5. 

see sort of a plateau in the correlation function in Fig. O 
which can be systematically decreased by increasing D. 
Our results indicate that the interaction does not destroy 
localization, at least for large disorder. This appears to 
be in contrast with the predicted metal-insulator transi- 
tion at a finite T in the corresponding fermionic system--- . 

However, in light of observed logarithmic growth of 
entanglement entropies with time we cannot exclude an- 
other possibility of very small diffusion constant and, in 
strict sense, absence of localization. Yet, another pos- 
sibility would be existence of localization delocalization 
transition at smaller disorder strength /i, but this regime 
is much harder to simulate with the present method as 
entanglement entropies increase faster for smaller h. 

Still, we were able to verify that the entropy grows 
logarithmically also for smaller fields, e.g. for = 2, as 
well as for other interaction strength, e. g. A — 1.5. For 



example, for A = 0.5 and h = 2 we find S{t) ~ clogt 
with c which is about 4 times larger than for h = 5, A — 
0.5. This indicates that the system probably displays the 
same kind of localization also for smaller magnetic field 
strength. If this is the case, transition in level spacing 
statistics observed i n^^i^^ might be just due to finite size 
effects (localization length being larger than the system 
size). 

We have also checked that the results of correlation 
functions, in the regimes which we show, do not signif- 
icantly depend on the system size n. For instance, the 
correlation function for n = 20 is practically the same as 
the one for n = 50 in Fig. [S] 



IV. CONCLUSIONS 

tDMRG is for general systems inefficient due to fast 
entanglement growth. Nevertheless, we have shown that 
it can be efficient for a disordered Heisenberg model, and 
potentially as well for some other disordered interacting 
models in one dimension. Spin-spin correlations evalu- 
ated in our study at large temperatures exhibit local- 
ization in spite of nontrivial interaction term, inferring 
that all-many body states are localized in large enough 
disorder. 

It should be noted that this is the first successful appli- 
cation of tDMRG for simulation of long time evolution at 
high temperature in a non-solvable (non-integrable) sys- 
tem. The efficiency of the method is related to an inter- 
esting observed logarithmic asymptotic growth of entan- 
glement entropy for time evolution of, both, pure states 
and operators. 
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